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Abstract 

We review the implementation of individual particle time-stepping for 
N-body dynamics. We present a class of integrators derived from second 
order Hamiltonian splitting. In contrast to the usual implementation of 
individual time-stepping, these integrators are momentum conserving and 
show excellent energy conservation in conjunction with a symmetrized time 
step criterion. We use an explicit but approximate formula for the time 
symmetrization that is compatible with the use of individual time steps. No 
iterative scheme is necessary. We implement these ideas in the HUAYNO0 
code and present tests of the integrators and show that the presented inte- 
gration schemes shows good energy conservation, with little or no systematic 
drift, while conserving momentum and angular momentum to machine pre- 
cision for long term integrations. 
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1. Introduction 

Astrophysical N-body problems often show a large dynamic range of 
timescales within their component systems. Instead of a fixed or varying 
global time step most current astrophysical codes update the force and ad- 
vance particles with a time ste p that is determined for each particle sep- 
arately (for a recent review see iDehnen and Readl . lioilh . Such individual 



particle time step schemes (IPTS) allow efficient simulations by concentrat- 
ing the computational resources on the parts of the system that experience 
the strongest evolution. The most popular IPTS is a block time step scheme 
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wher e the parti c les are organized in a pow e r of two hierarchy o f time step 
bins dMcMillanl . Il986l : iHernauist and Katd . Il989l : iMakinol . Il99lh . The bin 
in which a particle resides determines the frequency of force evaluations and 
each successively higher bin has a factor 2 smaller interval between updates. 
This block time step scheme has been adopte d in a wide range of codes , 
in fo r example codes for galacti c simulations (jDubinskil . Il996l : iMagorrianl . 



20071 ) . cosmol ogical simulations (IStadell. l200ll). Sm ooth Particle Hydrody- 
codes ( Springe! 20051 : IWadslev et al. . 2004) and also modern Her- 



namics 

mite integrators for collisional stellar systems ( Makino , 



Portedes Zwart et all 2001; lHarfst et all [20081 : iKonstantinidis and KokkotasL 



1991 



Aarscth, 1999; 



2010)) . The popularity of this scheme stems from the fact that it allows for 



individual tailored time steps while still grouping particles with similar time 
steps together - which means that the cost of synchronizing the rest of the 
system is shared by all the particles in a given bin and parallelization of the 
force calculation is possible. 

A notable exception where the simple block time step scheme is not 
used is when long term integration requires conservation of the integrals of 
motion to high precision. A prototype of a problem where this is required is 
the integration of planetary systems over billions of years. In this case care 
must be taken because IPTS typically do not conserve momentum and show 
drifts of the total energ y. For these problems symple ctic integrators based on 
Hamiltonian splitting ( Wisdom and Holman . 199ll ) are used, using the fact 
that the dynamics is dominated by a central object. For close encounters 
between objects orbiting the central body, additional splitting m ethods have 
been developed invq l ying the interaction terms betwe en orbiters (jChambersl . 



19991 : iDuncan et all Il998l : iMoore and Quillenl . l201lh . 



Here we will derive a class of IPTS integrators applying these ideas for 
general N-body problems by splitting the N-body Hamiltonian starting from 
an initial pivot time step. We subdivide particles in two groups, according 
to whether their time step is shorter or longer than the pivot. The time 
evolution of the slow particles is evaluated using the current time step, where 
the terms involving particles with a shorter time step are grouped into a new 
Hamiltonian to which we apply the same procedure twice with a halved pivot 
time step. This procedure is repeated recursively until the time evolution 
of all particles is evaluated. In this way different integrators can be derived 
(depending on the details of the split) that retain a close resemblance to 
the block time step scheme, but which will conserve integrals of motion, 
in particular the total (angular) momentum, to machine precision. We will 
show that, when used in conjunction with an approximately symmetric time 
step criterion, they show good energy conservation and do not show energy 
drift. In section[2]we present our derivation and our symmetrization scheme . 
In section [S]we present various tests of these integrators, while we summarize 
and discuss the results in section [H 
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2. Method 



2.1. Derivation 

The problem at hand is to evolve a set of particles X for a time dt under 
the dynamics generated by the Hamiltonian: 
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The Hamiltonian consists of momentum terms 

-,2 



Px^y^- 
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and potential terms 



Vxx 
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The evolution of the state of the system is given by the flow operator 
exp(dffl) where H is the Hamiltonian vector field corresponding to H. In 
cases the Hamiltonian can be split in two parts H = Ha + Hb, where 
expdtM-A and expdtMs can be calculated, the flow of H can be approxi- 
mated by composing expdtEU and exp riffle, e.g. 



exp(dtH) w exp(dt/2MB)exp(dffl A )exp(di/2H i j) 



(4) 



in th is case the approximation holds to second order in dt (jHairer et al 
200d ). An example of this is the division in Ha = Vxx and Hb = Px 



which generates the familiar second order Drift-Kick-Drift (DKD) verlet 
integrator. 

The derivation of approximate solvers by the division of the Hamilto- 
nian into simpler, solveable, subsystems is a general method not limited 
to the above division in momentum and potential terms. For example, 
efficient planetary integrators may be constructed by splitting the Hamilto- 
nian in a part that describes the Keplerian motion around the central star 
and a part that describes th e pertu rbation of the planetary bodies on each 
other ( Wisdom and Holman . 1 99 ll ) . Additionally, splitting can be used as 
a strategy for focusing the computational effort on the parts of the system 
that are evolvin g most strongly, either by decomposing the potential in two 
( Springe! 2005 . for a split into long range and short range for ces in cosmo- 
logical simulations) or more components (jDuncan et all Il998l . employing a 
series of succesive shells) , or b y grouping particles according to their dynam - 
ical evolution timescale (e.g. iFujii et al.l . 120071 ; ISaha and Tremaind . 11994 ). 
Here we derive multiple time step methods by choosing the division of the 
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Hamiltonian adaptively and recursively based on the current time step as- 
signed to the particles such that one system S contains the interaction terms 
of all the particles with a time step larger than dt, and the other system F 
contains all the interaction terms of the particles with a time step smaller 
than dt. This method of s ubdividing the system acco rding to the timestep 
is very close to splitting in ISaha and Tremainel (| 19941 ) (the "PASS" method 



we present below employs the same split as their method), apart from the 
fact that we allow the partitioning of the system to change continuously. 
This last property means that we have to be careful choosing the timesteps 
(discussion of which we defer to section 12, 2p . 
For example, we may choose the split 

H s = Ps + Vss + Vsf 

H F = P F + V FF (5) 

(Note that indeed H = Hs + Hp) and approximate 

exp(cMH) « exp(dt/2M F ) exp(dtM s ) exp(dt/2M F ) (6) 

The evolution operator exp dtMs is then approximated by the DKD (or any 
other second order scheme) . This consists of drifts of the particles in S (from 
Ps) and kicks on the particles in S and F (from Vgs + Vsf). These kicks 
consists of updates of the particle positions where every —Gniimj /\r{ — rj\ 
term in the potential part generates kicks 

dpi = dt- — |3 fa - rj) 

\ r i r j\ 

, Gmimj , 

d Pi = dt \r -r 3 ^-^) ( 7 ) 

For the two exp dt/2Mp operators of the F system we apply the same proce- 
dure as for the expdtW operator (with halved time step). The splitting can 
thus be applied recursively to the evolution operators with time step dt/2 k . 
The recursion ends when a time step sufficiently small is reached such that 
no particles end up in the F system (of that level). Note that we made a 
choice in subdividing the system: Eq. [5] is certainly not the only possible 
split. For example, the potential cross terms Vsf could also be put in Hp, 
resulting in a different method. The difference between the two is that in 
the former case (eq. [5]) the interactions between particles in S and F systems 
are done on the slowest of the time steps of the i and j particle pair. In the 
latter case the interactions propagate to the fast system, and will in the end 
be updated on the minimum of the time steps of each particle pair, 

thus the split in eq. [S] can be expected to be computationally faster. Note 
that the integration scheme that results is similar to conventional individual 
block time step schemes, where the time step is subdivided in powers of 
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hierarchy of time step bins to dictate the frequency in which the total forces 
for the particles in that bin are updated. A particle in a fast time step bin 
is evolved using a force calculated from all other particles - the particles 
that need less updates (higher in the time step hierarchy) are extrapolated 
to the time of the fast particle if necessary. In the scheme presented here 
the total force is never (or only incidentally) calculated: all the velocity 
changes are effected by momentum kicks (which need only originate from a 
subset of the particles) - and always pairwise: every i <— j interaction kick 
is accompanied by its corresponding i — > j interaction kick. The scheme 
is therefore manifestly momentum conserving. The extrapolating schemes, 
including their Hermite generalizations, are not momentum conserving. 

By putting the Vsf m the slow system in Eq. [5] we made the choice 
of calculating the interactions between 'slow' and 'fast' system on the time 
step of the slow particles - hence we refer to this integrator as the HOLD 
integrator (the Vsf axe 'held' on the slow time step). Alternatively, when 
we put the Vsf terms in the fast system: 

H s = Ps + Vss + Vse 

H F = Pf + Vff + Vfs + Vfe (8) 

Note that we have included the Vse and Vfe-, these denote any interactions 
that are inherited from the higher levels. As we will see below this integrator 
does not conserve the position of the center of mass, the corresponding 
integrator that does is 

Hs = Vss + Vse 

H f = Ps + Pf + Vff + Vfs + Vfe (9) 

We refer to this scheme as PASS (and to eq. [8] as PASS1). 

We will compare the HOLD and PASS integrators with the conventional 
integrator without individual time steps but with shared adaptive time steps 
(SHARED) and the conventional block time step integrator which extrap- 
olates particle positions (BLOCK). A cartoon representation and overview 
of the differences between the different integrators is given in Figure [TJ 

2.2. Symmetric time- stepping 

The above integrators are symplectic only as long as the particles do not 
change their timestep. For general N-body simulation codes, where parti- 
cles experience very different dynamical timescales during their evolution, 
we need to relax this condition and let the particles move between timestep 
sets. Using adaptive timesteps dependent on the phase space coordinates 
will in general destroy th e symplecticity ( and hence the conservation prop- 
erties) of the integrator ( Skeel and Gearl . 19921 : Preto and Tremaine . 1999 : 
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Figure 1: Schematic overview of the different integrators. Panel A shows the interactions 
that are calculated between two particles (1 and 2) by the SHARED integrator for a 
fiducial trajectory of these two particles, where particle 1 has a time step 4x larger than 
particle 2, which is also globally in the smallest time step bin (with the rest of the system 
is represented by the greyed out trajectory). All the interactions between the particle are 
calculated on the smallest time step. Panel B shows the case of the conventional block 
time step scheme (BLOCK integrator). In this case the total force on each particle is 
calculated determined by the individual time step of that particle. Panel C shows the 
HOLD integrator (Eq. [5]). In this case particle 2 only feels the interaction of particle 1 
with a frequency determined by the time step of particle 1. In between the particle may 
be kicked by interactions with other particles. Panel D shows the PASS integrator (Eq.[9j). 
For this integrator the mutual interactions between 1 and 2 are calculated on a time step 
determined by the fast particle 2 (note that in this figure the same number of interactions 
is calculated for the SHARED and PASS integrators - in general this is not the case). The 
difference with panel B is the fact that each kick of the velocity of particle 2 is matched 
by the corresponding kick on particle 1. 
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Mikkola and Tanikawal . [l999), although with some restrictions it is possi- 
ble to construct adaptive timestep symplectic integrators by cons i dering 

the Hamiltonian in an exten ded phase space (jPreto and Tremainel. Il999l: 

Mikk ola and Tanikawa , 19991 ) or using a variational approach (jFarr and Bertschinger 
20071 ). For the integrators we present here it is however possible to recover 
long term energy conservation by ensuring time reversibility of the integra- 
tor. This can be done by using a time-symmetrized time step criterion for 



the particles ( Hut et al. . 1995 : Preto and Tremaine . 19991 ) . 



To be specific, we start by reviewing in more detail why the numerical 
integration of a particle trajectory loses its time symmetry (and hence its 
energy conserving properties) when a variable time step is adopted. Given a 
particle trajectory (r(t),v(t)) one tries numerically to integrate this using a 
time step function r(t) (which is a function of time through its dependence 
on the state variables r and v). If we advance a state r(t) — > r(t + T\), v(t) — > 
v(t + Tx) we reach a new state with a new time step T2- If we would run the 
integration backwards in time, in general t\ ^ T2, so we would not arrive at 
the original state, and due to the numerical errors we would not follow the 
same numerical approximation to the trajectory. This situation would not 
arise if the time step function had the property that r(i + r(i)) = r{t). Note 
that r(t) refers to a time step in positive direction while r(i + r(t)) refers to 
a time step in the negative direction, so the requirement should be read as 
t — (f + r(t)) = r + {t). In fact we can introduce such a symmetric time step 
function r sym : 



VW = '■(t)/2 + r(i-r^ n )/2 

r+ m (t) = r(t)/2+r(t + r+ m )/2 (10) 



this becomes a workable ansatz if we approximate 



mi 



so we can write: 



r(t) 



i sym 



V 1 2 dt > 



(12) 



which is the approximation we will use below. For a time step proportional 
to the inter-particle free-fall times 




(13) 



with mj = rrii + rrij , the derivative is 



dt 



3V{j ' T{j 

2rl 



(14) 
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hence the symmetrized version of the free-fall time step criterion becomes 



" " ■■ ' (15) 



2 dt I J 

A time step proportional to the interparticle flyby times 



Tij=m— (16) 



with _ _ 

dTij • Vij . jiij 



dt 



2 v \ 1 2 
r ij v ifij 



1 + (17) 



can be similarly made symmetric (note that eq. [T7| simplifies the expression 
by only considering the i and j particle for the time derivative of Vij). We will 
take the particle time step to be equal to the minimum of the symmetrized 
free-fall (defined by eqs. Q21 [EH and and symmetrized flyby time step 

(eq. ma in and Hsu . 

T here are various other ways in the literature to symmetrize the time 
step. iHut et al.1 (jl995h advocated an iterative procedure to symmetrize Eq. 



[TU1 In this case an trial step is taken, after which the average of the new and 
old time step is taken. This can be repeated to converge to the solution of 
Eq. [TOl This has the drawback that the iteration is an expensive operation 
which needs multiple force evaluations for the system. Alternatively, a time 
step determined by 

ToldTnew = r(t) 2 (18) 

provides an explicit symmetrized time step (jDehnen and Read! 1201 ill. Both 



of these do not work when used with a block time step scheme (jMakino et al. 



2006 : Dehnen and Readl . 2011 ). For the former iterative solution this is 



because of a "flip-flop" problem where convergence is not reached because 
of the restriction of the time step to power s of tw o of a base time step. 
A solution in case is available ( Makino et al. . 20061 ). but entails saving the 



whole system state over a time interval and multiple iterations. 

To test whether the expressions given here (which due to the truncation 
after first order in r in Eq.[12]are only approximately time-symmetric) yield a 
conservative variable time step method when using a power of two hierarchy 
of time steps, we integrate a Kepler orbit with eccentricity e = .9 for 1000 



orbit s and a ry = 0.01 (this is the same test as fig. 1 and 2 of IHut et al 



19951 ). The test uses the symmetrized free-fall time step (eq. [T5|) . The time 



steps are chosen to be block time steps as in the integrators presented in 
section 12.11 (for the two body problem the integrators of section I2TT1 and the 
usual block time-stepping are equivalent). The orbital elements for the orbit 
are plotted in Fig. [5] . As can be seen there is no systematic drift in energy 
(the orbit does precess with a constant speed). The performance of the 
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integrator is similar to the iterative scheme of (jHut et all Il995l ) , and does 
not show the systematic drift in a and e of a non-symmetrized integrator. 

In Figure [3] we show the resulting (maximum) energy error for a set 
of 10 4 orbit integrations for different values of the eccentricity and time 
step parameter 77, plotted vs the number of steps per orbit. Although the 
integration scheme is in general robust, this figure shows that for a given 77 
it tends to keep the number of steps per orbit constant. This could be a 
disadvantage if such a time step is used in a run with multiple bodies where 
encounters with different e can occur. This can be remedied by choosing an 
error control time step r ec 

Tec = r(-y (19) 
TO 

with 7 a small exponent such that the time steps become smaller as r de- 
creases. Such an expression is also easily symmetrized. Figure [3] also shows 
the results for this time-stepping scheme (with 7=1/3). As can be seen, in 
this case runs with constant 77 show a much less marked increase in error, 
showing an almost flat profile for the orbits with extreme e (albeit at the 
price of a rapidly increasing number of steps per orbit) For the other tests 
below we do not use the r ec as these tests are not especially sensitive to the 
high eccentricity behavior. 

2.3. Implementation 

The integrators are implemented in the HUAYNO (for Hierarchically 
split-Up AstrophYsical N-body sOlver) N-body code, released as part of the 



Astrophysical Multi-purpose Software Envrionment (AMUSE JPortegies Zwart et al 



2009 )q HUAYNO is written in C as testbed for the N-body integrators de- 
scribed here. There is a very close correspondence of procedures in HUAYNO 
and the mathematical operators for the algorithm: evolve, drift and kick 
operators have their counterpart in the code (this is illustrated in the pseu- 
docode for the HOLD integrator given in fig. The particles themselves 
are stored in an array which is reshuffled to match the partitioning required 
by the integrator (this partitioning is very sim ilar to the block time step 
implementation of McMillan and Aarseth . 19931 ). The subsystems are then 



contiguous memory blocks conveniently referred to by a pointer and size. 
This has the added benefit that the particles which end up accessed the 
most end up close to each other in memory. The shuffling itself is an 0{N) 
operation with negligible impact on the performance. 

The code implements an optional Plummer softening parameter e for the 
gravitational interactions. It allows to set separately 771 and 772, for the time 
step parameter of the free-fall (eq. [13|) and flyby (eq. [T6]) time step criteria. 
If not explicitly mentioned they are taken to be the same (and referred to by 
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orbit 



Figure 2: Orbital elements for the variable time step leapfrog integrator with symmetrized 
time step criterion feqs. 1131 1141 and I15[) . Plotted are the changes in semi-major axis a, 
eccentricity e, longitude of apo-center 1 and the time of apocenter passage as a function 
of orbit number (sampled at apocenter). The time steps are varied in the same way as in 
the blocktime step schemes. 
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Figure 3: Maximum relative energy error vs number of steps per orbit. For this a Kepler 
orbit is integrated for 10000 orbits. The lines in black use the standard symmetrized time 
step criterion (eqs. 1131 1141 and !15p for different values of the eccentricity (from bottom to 
top e = 0.6,0.9,0.99,0.999,0.9999,0.99999), dotted lines connect runs with the same r). 
Thin red lines show the same for runs with the error control time step Eq. 1191 



Figure 4: Pseudo code of the HOLD integrator, 
function evolve_hold(system, dt) 
calculate_timestep( system) 

fast , slow=sort_in_f ast_and_slow(system,dt) 

evolve_hold(f ast , dt/2) 

drift (slow, dt/2) 

kick(slow, fast, dt) 

kick(fast, slow, dt) 

drift (slow, dt/2) 

evolve_hold(f ast , dt/2) 
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Figure 5: Error in conserved quantities for N = 100 plummer sphere runs. Shown is the 
energy error (upper left panel), drift in the center of mass position (upper right), linear 
momentum (lower left) and angular momentum (lower right) for different integrators: 
SHARED (thick black line), BLOCK (thin red line), PASS1 (dashed thin red line), PASS 
(black dashed line) and HOLD integrator (dotted black lines). 

77) in the tests below. Unless explicitly given the tests below use unsoftened 
gravity, e = 0. 



3. Tests 

To examine the differences between the 5 integrators presented in section 
12.11 we evolve an N = 100 Plummer sphere of particles with equal mass. 
Total mass is 1 and uni ts are scaled such that G= l, total mass M=l and 
E=-0.25 (Nbody units, iHeggie and Mathieul . Il986h . Figure [5] shows the 



evolution of the Energy E, momentum P, angular momentum L and the 
center of mass position X for each integrator, all evolved with the same 
77 = 0.01. 

The first thing to notice in these plots is the different behavior of the 
BLOCK integrator with respect to momentum conservation. The error in 
center of mass and momenta (linear and angular) remain at machine preci- 
sion for the conservative integrators, whereas the BLOCK integrator shows 
an error of ~ 10~ 7 . Note also that the center of mass of the PASS1 in- 
tegrator run is not conserved - this may be suprising: this integrator does 
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Figure 6: Growth of relative energy error for 20 symmetrized (black lines) and 20 non- 
symmetrized time step (red lines) runs. This is the same test as in Makino et al. (2006). 
The same set of N = 100 Plummer sphere initial conditions were used for the two sets 
of simulations. Symmetrized runs took w 20% longer than unsymmetrized runs (26 vs 22 
sec in total). 



conserve linear momentum because of the pairwise force kicks (as can be 
seen in fig. [5]) . The reason for the shift in the center of mass is the fact that 
this integrator includes kicks between particles that are not synchronized 
in time: combining Eq. [6] with Eq. [8] we see that midway in the half step 
of the evolution of exp dt /2Hp the kicks between the S and F system are 
calculated. At that point the positions of the S system particles lag dt/A 
with respect to the F system. After the kicks are applied the position of an 
S particles drifts an extra Avidt/4 from the first Av± kick. In the second 
half step this is compensated by missing a Av2At/4 drift, but these two will 
in general not cancel out exactly: hence the center of mass position is not 
conserved. This will happen irrespective of whether a DKD or KDK split 
is taken as the basis for PASS1. The proper formulation (PASS), which en- 
sures time synchronized kicks, does show conservation of the center of mass 
position (fig. [S]). Note also that applying the KDK method with the HOLD 
integrator (for exp dtWs in eq. [5]) will result in kicks between particles at 
different evolve times, and thus a deviation of the center of mass. 

The energy errors are very similar for all the integrators. They remain 
all at a level of fewxl0~ 7 and do not show any linear drift. The reason that 
the BLOCK integrator performs competitive with the other integrators is 
that the errors in momentum etc. remain small enough not to induce energy 
errors (Note that the Plummer sphere initial conditions represent a problem 
that is not especially sensitive to errors in momentum conservation). 

As we mentioned above, our symmetrization scheme (eq. I12|) of the time 
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timestep parameter time(s) 

Figure 7: Energy error as versus time step parameter and run time. Each run evolves an 
N = 1024 plummer sphere for 1 N-body time unit (lines use the same legend as in figEJ. 
Left panel shows relative energy error as a function of time step parameter 7/, right panel 
shows error versus the wallclock time elapsed for that run. Energy error is the maximum 
recorded over the run. 



step is not exact, approximating the formal symmetry of Eq.llOl The Kepler 
test of the symmetrization scheme in section [2T2l suggests that it is adequate, 
however for larger N simulations the additional terms in the timestep deriva- 
tives in EgfTZl or the selection of minimum timesteps in Eq. [15] could cause 
the timestep symmetrization to fail. Hence, we conduct a n additional test of 
our scheme by repeating the test of lMakino et all fxxxt ) presented in their 
figure 4. This test consists of running N = 100 Plummer sphere m odels for 
symmetrized and non-symmetrized time steps. Makino et al. ( 20061 ) showed 
that the linear drifts of the non-symmetrized time step runs are suppressed 
when using the symmetrized version (they used a scheme for block time steps 
using 6 iterations to symmetrize the time step). Our result for this test is 
shown in figure The integrator we use for b oth sets of run s is the HOLD 
integrator. The result is very similar to the Makino et al. ( 20061 ) figure. 



This shows that our symmetrization scheme works at least as well as the 
approach given there, and seems sufficient for the second order integrators 
employed here. 

In figure [7] we show the behavior of the integrators as a function of the 
time step parameter rj. All integrators show the second order dependence 
on rj as expected. The actual energy errors at the same rj also lie within 
a narrow band. If we look at the trend with run-time we see that the 
SHARED integrator is a factor 300-1000 slower than the other integrators. 
Within the integrators with individual time steps the HOLD integrator is 
most efficient, then the PASS1 and PASS integrators and lastly the BLOCK 
scheme. Looking at energy error dependence on rj we see that BLOCK 
and HOLD are actually very similar, so the difference (of a factor ~ 3) in 
efficiency as measured by wall-clock time is a result of the lower number of 
force evaluations of the HOLD method. 
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Figure 8: Lagrangian radii, energy and momentum momentum error and relative run- 
time for the N=1024 Plummer sphere core collapse test with smoothing. Upper left panel 
shows Lagrangian radii for runs with the HOLD (dotted black lines), PASS (dashed black 
lines) and the BLOCK (drawn red lines) integrators (lines in the other panels use the 
same coloring). Upper right panel shows the relative energy error, lower left panel shows 
the momentum error and lower right shows the relative cpu time per unit step for the 
simulation with respect to the HOLD run. 



3.1. Core collapse 

As a last test we evolve a N = 1024 Plummer sphere of equal mass 
particles to core collapse. First we employ smoothing with an e = 1/256 
(r/ = 0.01). In this case we evolve to t = 700. The setup was chosen to match 
the test in Nitadori and Makino ( 20081 ). We plot the Lagrangian radii, en- 
ergy and momentum error in figure [5] (The angular momentum and center 
of mass errors are omitted as they show the same pattern as the momentum 
error). We plot the results for the PASS, BLOCK and HOLD integrators. 
As can be seen in the plot of the Lagrangian radii all integrators show the 
expected evolution of the mass distribution. Differ ences within the runs 
fall w ithin the expected statistical variation (see also iNitadori and Makind . 
20081 ). Note that although the runs start with the exact same initial con- 



ditions these can be considered effectively different runs as they diverge 
exponentially on the crossing time timescale. The energy errors in Figure [8] 
show similar results for the integrators, the BLOCK integrator showing a 
consistently factor 2 to 3 larger error. Note the absence of energy drift also 
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Figure 9: Lagrangian radii, energy error and wallclock time for the N=1024 Plummer 
sphere core collapse test with e = for the HOLD integrator. Upper panel shows la- 
grangian radii, Middle panel shows the relative energy error, lower panel shows the wall- 
clock time for the simulation. 
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in this long term integration (the energy errors do increase somewhat as the 
integration progresses, but they return to zero and change sign). In the plot 
of the momentum error we can see that the PASS and HOLD integrators 
conserve momentum to machine precision. A slight increase in error is visi- 
ble at the end. This is due to the loss of precision in calculations involving 
escaping stars. If we compare the time needed to advance the simulation 
of the three integrators we notice a clear difference: at the beginning the 
HOLD integrator is a factor ~ 3 faster (which is consistent with the tests 
in fig. [5]), as the simulation evolves the difference between HOLD and the 
other two integrators increases to a factor 12 at the end. The HOLD inte- 
grator becomes more efficient as the timescale differences within the system 
become larger. 

As an additional test for the HOLD integrator we run the same prob- 
lem for e = 0. This is a demandin g test as in this c a se str ongly bound 
binaries form during core collapse (jGiersz and Heggid . 1 19941 ). This test 
calculation is shown in figure [9l In this case the core collapse occurs at 
t ~ 350, and the evolutio n of the Lagrangian radi i and core radius shows 
a pattern consistent with Giersz and Heggie ( 19941 ). The energy errors in- 
crease during the core c ollapse (for reference this can be c ompared with 
current Hermite codes in Konstantinidis and Kokkotasl . 2010 . their fig. 21). 
Note that the error increases during very short episodes, after which the 
error remains bounded. These error jumps occur during strong interac- 
tions with the binary that forms in the center, and may be caused by 
failures of the timestep symmetrization, although similar error jumps oc- 
cur also in codes using compl e tely d ifferent timestepping criteria (see again 



Konstantinidis and Kokkotasl . [201CJ, for an example) A more uniform error 



behavior may be obtained by decreasing the r/ parameter during these 'rough 
patches,' for a modest increase in computation time (for very close encoun- 
ters the loss of accuracy may be due to the numerical precision, in this case 
some form of regularization is necessary). Plotting the wall-clock time of 
the simulation we see that towards the end of the simulations slows down 
(for the BLOCK and PASS integrators this would mean an additional slow- 
down by a factor 100). This is due to the formation of hard binaries. The 
HOLD integrator therefore is able to integrate systems without smoothing, 
although a more efficient handling of binaries is desirable. 



4. Discussion 

We have presented a set of conservative individual time step integra- 
tors derived using recursive Hamiltonian splitting. These integrators are 
remarkably simple and closely follow previous work: specifically they can 
be seen as a generalization of the scheme used in the mes h-based integra- 
tor for collision-less systems GROMMET ( Magorrianl . 120071 1 and are closely 
related to the Hamiltonian splitting methods used in planetary dynamics 
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([Duncan et all Il998l : ISaha and Tremainel. Il994f) or the force s plitting meth- 



ods used in molecular dynamics (e.g. Tucker man et al. . 19921 ). In contrast 



to the usual formulation of block time steps they conserve momentum and 
angular momentum to machine precision and in conjunction with the (ap- 
proximately) symmetrized timstep criterion presented here they show good 
conservation of energy. Our tests show that the energy conservation is a con- 
sequence of the use of a symmetrized time step criterion (as the conventional 
BLOCK integrator shows similar error behavior with the symmetrized time 
step). Of the new integrators HOLD and PASS the former is more efficient 
because it calculates the interaction between two particles with different 
time steps at an interval determined by the slower of the two particles. The 
Hamiltonian splitting considered here can be applied more generally. For 
example, Smooth Particle Hydrodynamics (SPH) codes also use a similar 
time-stepping scheme, and the scheme presented here directly carries over 
to SPH codes. 

One drawback of the subdivision in terms of the time steps assigned to 
the particles is that the interactions for a given particle in a time bin are 
calculated together with all other particles with the same time step. For 
example, if the system contains multiple binary stars (and assume all these 
binaries have the same orbital parameters) then the stars that form these 
binaries all have the same time step and will end up in the same bin, with 
all their mutual interactions calculated on the fastest time step. The worst 
case scenario (which is however physically significant) is that all the stars 
are in binaries. In this case the split time step scheme as presented here does 
not save computing time with respect to a shared time step simulation (this 
is also the case for conventional block time step integrators). To ameliorate 
this one would need a subdivision at each level of the system in multiple 
systems based on e.g. a spatial partitioning at each level, as opposed to a 
simple split in a fast and a slow system. 

We have limited ourselves to second order integrators. The second or- 
der nature of the integrators is a consequence of Eq. this is the second 
order leapfrog composition. The re is no reason w hy we could not use a 
higher order composition scheme ([Mclachlanl . Il995l ) , as long as we take care 
selecting the same (or at least of the same order) composition schemes for 
the component integrators. A problem with using higher order integrators 
is that the compositions contain terms with negative time steps. While in 
principle this is not a problem, it does present a problem for the efficiency of 
the integrator as the time interval that the integrator follows is bigger than 
the integration time interval by a factor a > 1. Note that within a recursive 
scheme every level in the time step hierarchy will increase the length of this 
interval by this factor, so that unless a is very close to 1, a k can increase 
rapidly. 

In addition to the algorithmic enhancements mentioned above, the nu- 
merical implementation of the integrators can be improved. For large scale 
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production work they would need to be parallelized and probably also be 
adapted to use Graphic Processor Units (GPUs). At the moment a partial 
implementation on GPUs is available, where only the kicks are implemented 
on the GPU. This implementation suffers from the inefficiencies of transfer- 
ring data to and from the GPU. With these improvements we expect that 
the code will be competitive with modern Hermite codes. 
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